JOURNAL OF RESEARCH of the National Bureau of Standards— B. Mathematics and Mathematical Physics 

Vol. 71 B, No. 4, October-December 1 967 

Solving Equations Exactly 

Morris Newman 
Institute for Basic Standards, National Bureau of Standards, Washington, D.C. 20234 

(August 22, 1967) 

A congruential method for finding the exact solution of a system of linear equations with integral 
coefficients is described, and complete details of the program are given. Typical numerical results 
obtained with an existing program are given as well. 
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1. Introduction 

The problem of the solution of a given set of linear 
equations Ax = b on a high-speed digital computer 
has been studied intensively, and there are a large 
number of methods, more or less satisfactory, for 
carrying out such a solution. Nevertheless occasions 
arise when existing methods are inadequate, either 
because the solutions are required exactly, or because 
the coefficient matrix A is "ill-conditioned." A no- 
torious example of the latter is furnished by the Hilbert 
matrices A — H n given by 



H n 



1 



i+y-i 



i,j- 



Here even for moderate values of n (say 7 or 8) the 
gokitiog by any of the usual mclliuus becomes awk- 
ward, if not impossible. For these reasons, as well as 
many others, the method of solution described in this 
note is of interest. It is not at all sensitive to the con- 
dition of A, since it determines the exact solution (and 
not an approximate one) by number-theoretical 
methods. Thus the usual ills caused by round-off, 
truncation, etc., do not exist. It can fail, of course, 
but not for any of the reasons which cause the usual 
methods to fail. Its principal disadvantages are that 
it is limited to systems with integral elements, and that 
it is somewhat time-consuming. It is not particularly 
suited to hand computation, and definitely finds its 
role in high-speed digital computation. 

A related discussion and elaboration are given by 
I. Borosh and A. S. Fraenkel in their article. 1 



1 I. Borosh and A. S. Fraenkel. Exact solutions of linear equations with rational coefficients 
by congruence techniques, Math. Comp. 20, 107-112 (1966). 



2. Description of the Method 

Let A be a nonsingular integral n X n matrix, b an 
integral n X 1 vector. Put 

d=det (A) 

and denote the adjoint of A by A adj , so that A (l(lj is also 
a nonsingular integral n X n matrix, satisfying 

AA* d *=A & UA = dI. 

For any matrix B= (bjj) define 

M(B) = max \by\. 
i.j 

Let x be the solution of the system 

Ax=b. 

Then 






1 
d 



A^b = -jy, 



where 



y = A^b 



is also an integral vector. 

The method is based on the following simple obser- 
vation: Suppose that m is an integer such that 

(d, m) = l, 

m>2 max(|d|, M{y)). 

Then any pair of solutions d m , ym, of the congruences 

d = d m mod m, 

Ay m = db mod m 
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satisfying 



\d m \ < hm, 
M(y m ) < \m 



must in fact coincide with d, y. 
It is clear that d= d m , since 



\d\ < 2/n, \d m \ < 2 m. Furthermore 



d = d m mod m and 
since Ay = db 
mod tti, Ay m — db mod m, we have A(y — y m )—0 
mod m. Hence A a<ii A(y — y m ) = mod tti, d(y — y m ) 
= mod tti; and since (d, 77i) = l, it follows that 
y = y m mod tti. Finally, since both M(y) < 2m, 
A/(ym) < 2 tti, we have y = y m . 

In applying the observation above we choose 
m = m 1 • m 2 . . . m s , where (mi, m,j) = l for i#jf, and 
use the Chinese remainder theorem. Define ml by 



— ml = 1 mod /ft*, ^ m.' < m*, 1 = / = 5. 

TTlj 

Then the solution of the system 

-z = at mod mi, I t=z i = s 

is given by 



z = > — 711 . at mod 
• , mi 



For each i, I ^ i 



determine d n 



l V ym i 



that 



d= d m . mod TTij, 

Ay m . = db mod m,. 

(We shall show how this can be done without having 
computed d previously.) 
Next determine d m , y m by 



d m — V — m'. d m . mod tti, 
& m l l 



y m = V — Tri^y™. mod m. 

Finally, determine d, y so that 

d = d m mod m, 
y= y m mod m, 

and |rf| < 277i, M(y) < im. 

The practical utility of this method rests primarily 
on the fact that the bulk of the computations are per- 
formed modulo the "single-length" numbers ttij, 
1 ^ i ^ s. The final reconstruction according to the 
Chinese remainder theorem is done accumulatively 



throughout the computation and is the only time when 
multi-length operations are required. In practice the 
moduli mi are chosen as prime powers, for then the 
solution of the congruential system 

d= d m . mod Tiii, 

Ay m = db mod 771 j 

becomes particularly simple. 

The ideal program would first determine a permis- 
sible value for tti. This can be done by Hadamard's 
inequality, for example, which states that the absolute 
value of a determinant does not exceed the product 
of the euclidean lengths of its row vectors. This im- 
plies that 

\d\^J M(A)\ 

M(y) Sn(*-lp~ M(A) n -m{b) 

which are at once derivable from Hadamard's ine- 
quality and from the specific form of the elements of 
A adi as cofactors of the elements of A. Thus it is cer- 
tainly sufficient to choose 

n n-l 

m >c= 2 max (n*M(A) n , n(n-l) 2 M(A) n ~ l M(b)) 



and satisfying 



(m, d) = 1. 



The simplest way to do this effectively would be to 
generate a sequence of different primes pi, p 2 , . • -,Ps 
such that (d, pi) = 1, 1 ^ 1 ; 2| s, pip 2 ■ • . Ps> c; and 
then choose m = pip 2 . . . p s . 

In practice, the procedure outlined above is unneces- 
sarily conservative and time consuming. A good prac- 
tical alternative is to operate instead with a prede- 
termined set of moduli ttij, knowing full well that in 
certain instances the process will fail. Failure prin- 
cipally occurs if (mi, d) > 1 for some i. By choosing 
the TTij as large primes, the probability of such an oc- 
currence can be made so small that this is not an 
important practical consideration. 

3. The Detailed Program 

Suppose that the machine for which the program is 
being written is capable of multiplying or adding 
together any pair of whole numbers < K in absolute 
value. We choose as moduli the s largest primes 
Til* < K, and precompute and store the numbers 

m 

— , ml. Notice that ^ tti- < mi < K, but that the num- 

mi l l 

m 
bers — must be stored as multilength numbers. In ad- 

TTli & 



dition in and 



are required and must be stored 
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similarly. The choice of 5 is a matter of expediency: 
If K ~ 10 10 then 5=10 should suffice for most pur- 
poses. The program (in broad outline) is as follows: 

(1) Read in n, A, b. 

(2) Clear s-length answer cells d, y. 

(3) Set i=l (modulus tally). 

(4) Transfer A, b modulo m, to working temporaries. 

(5) Determine d m ., y m . so that 

ii 

d m . — d mod rrii, 
Ay m . = db mod mu 

and d m . and the components of y m . lie between and 
mi— 1, inclusive. (As will be shown, d need not be 
known.) 

(6) Accumulate s-length results: 



d-\ m'idnu ~* d, 



in , 

y+—m i ym i -*y. 

mi l l 

(7) If i < s, replace i by i+ 1 and go to (4); otherwise 
go on. 

(8) Reduce d, y modulo m so that 

\d\ < 2 m, 
M(y) <h m. 

(9) Check 5-length answers by multiplication: 
Ay=db. 

(10) Compute (in floating) the components of 
1 

X = d y ' 

(11) Print results (d, y as 5-length integers; x as 
floated numbers). 

The basic computation of course takes place in 
(5), and we describe this in some detail. Suppose then 
that p is a prime, and we are interesteH in solving 

d p = d mod p, 

Ay p = db mod p. 

Let B=(bij) be the nX(^+l) matrix (A, b) . The 
procedure is as follows: 

(5.1) Set k= 1. 

(5.2) Set S= 1 (cell in which d p is calculated). 

(5.3) Determine r so that r ^ kfkn, and {b r k, p) = 1. 
(If no such r exists, the machine prints out the infor- 
mation that the system is singular modulo p and halts.) 

(5.4) If r= A\ go to (5.6); otherwise go on. 

(5.5) Interchange rows r and k, and replace S by 
— S. Continue to denote the resulting matrix by B. 

(5.6) Determine b' kk (by the euclidean algorithm) so 
bkkb' kk = 1 mod p. 



(5.7) Replace S by bkkS, and bkt by b' k J>ku kt^f^kn 
+ 1, all computations being performed modulo p. 
Continue to denote the resulting matrix by B. 

(5.8) For 1^5^m, s ¥" k, replace b s t by 



bskbk 



k^t^n+1, 



all computations being performed modulo p; and con- 
tinue to denote the resulting matrix by B. 

(5.9) If k < n, replace k by /c+1 and go to (5.3); 
otherwise go on. 

(5.10) d p = S. 

(5.11) vf - Sb in+ i (mod p), l^i^n. 

It will be seen that the time required to accomplish 
step (5) is approximately equivalent to that required 
to solve a single linear system by the elimination 
method. Thus the total time is roughly that required 
to solve 5 such systems, and so is quite moderate. 

We note one or two other points: In the accumula- 
tion described by (6), we reduce m'.d m . and m'.y m . 

m 
modulo m\ before multiplication by — . The effect of 

this is that at the end of step (6), ^ d < sm, M(y) < sm. 
The reduction required by step (8) is then relatively 
easy to accomplish. 

It is certainly possible to provide for the contin- 
gency that (d, mi) > 1 for some i by a more elaborate 
program. If this occurs, the faulty m\ must be discarded 
and a new one substituted. The necessary constants 
would then have to be recomputed and the problem 
started anew. 



4. An Existing Program and Numerical Results 

A pilot version of this program was written for the 
Q32 computer at the System Development Corpora- 
tion, Santa Monica, Calif. The computer was used on 
a "time sharing" basis, under partial support of the 
National Institutes of Health on project number 
2050404. The writer would like to thank Russell Kirsch 
of the National Bureau of Standards for generously 
providing access to this computer. 

The following 10 primes were chosen as moduli: 
9999889, 9999901, 9999907, 9999929, 9999931, 9999937, 
9999943, 9999971, 9999973, 9999991. 

In the appended JOVIAL program, the lines have 
the following function: 

1-8. Array declarations. 
9-26. Precomputed constants. 
27-40. Read in n, A, b. 
41-42. Clear answer cells. 
43-44. Set up modulus p. 

45-51. Transfer A, b modulo p to temporaries. 
52-74. Solve Ax p = b mod p. 
75-78. Compute d p , y p . 
79-85. Accumulate d, y. 

86. Reset modulus. 
87-93. Reduce to range (—2 m, 2 m). 
94r-97. Print exact results (y, d). 
98-109. Multilength normalization subroutine. 
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110-117. Euclidean algorithm. 

118. Error print. 
119-122. Compute and print floated results (x). 



A A\ A 2 A 3 A 4 

A 5 A« A 7 A 8 A 9 



A number of test problems were run with uniform stand for the 10-length integer 
success, among which were the Hilbert matrices of all 
orders up to and including n=13 (a limit imposed by 
considerations of time). Define 



* n = L.C.M. (1,2,3, . . ., 2n-l). 

Then t n H n has integral elements. We chose b as the 
first unit vector, so that the output was the first column 

of {t n H n )~ l = — tin 1 . Since Z/,7 1 is known exactly, 2 we 

were able to verify the results. In the numerical results 
that follow, the lines 



^o+lO^-h. . . -hl0 6 M 9 . 
The print-out has the form 



y 
d 



We give the results for n = 9, 10, 11, 12 as repre- 
sentative examples. 



2296512 8395179 399 



-1860480 -5807169 -15993 



8876160 9525337 205250 



-3256960 -7152027 -1231505 



8085120 5744088 4002393 



-1092224 -6722299 -7471134 



6170240 1488177 8004787 



-2097280 -850387 -4574164 



3134080 9645230 1080010 



5508480 132897 480614 6 



+.661103602E-005 

-264441441E-003 

+.339366516E-002 

-203619910E-001 

+.661764706E-001 

-123529412E + 000 

+.132352941E + 000 

-756302521E-001 

-K178571429E-001 

1454720 2473993 2532211 4345 



-7008640 -7462660 -344456 -215090 



2138240 9402571 5511307 3441440 



-7048320 -8802831 -8460751 -6097590 -2 



1602944 8971893 5535157 9609881 10 



n = 9 

* w = 12252240 



*r!h ^- Savage and E. Lukacs, Tables of inverses of finite segments of the Hilbert matrix, 
NBS AMS 39, 105-108 (1954). 
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n = 10 

-9007360 -7429732 -8837894 -4024703 -27 

*„ = 232792560 

2773120 845307 5372030 7561453 41 



-1275520 -3984119 -8999527 -279180 -38 



5637760 6992059 4499763 139590 19 



-2412416 -2253879 -2061061 -140580 -4 



7928832 7977568 9387428 5426211 1011 



+.429566993E-006 

-212635662E-004 

+.340217058E-003 

-257997936E-002 

-K108359133E-001 

-270897833E-001 

-K412796698E-001 

-375939850E-001 

+.187969925E-001 

-396825397E-002 

7570112 599090 7074703 1 



-4206720 -5945445 -4482183 -102 



7031040 936185 7402580 1997 



-5623040 -2071066 -5757414 -18645 



4520960 5873099 2726424 97889 



-2467072 -6793918 -6724558 -313245 



1183360 416263 5752083 633949 

n=ll 

-8664320 -1963766 -252678 -815078 *„= 232792560 



1025920 2387982 1033370 645270 



-6011520 -3283547 -7125942 -286786 



6329472 7172313 1905861 54750 



9768320 338817 4351703 3285011 



-K519776062E-006 

-311865637E-004 

+.608137992E-003 

-567595459E-002 

+.297987616E-001 

-953560372E-001 

+.192982456E + 000 

-248120301E + 000 

+.196428571E + 000 

-873015873E-001 

+.166666667E-001 

126080 3885341 3441107 9248615 3937 



-4014720 -2801882 -1039178 -6275997 -281561 
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3676800 8710589 914159 3106599 6569771 



-3864000 -5494130 -7784298 -2449239 -3909927 

-70000 

4729600 7162434 5819510 3675134 3023534 

47 

-4094080 -1744323 -5084059 -3244700 -6326686 

- 187 

6241920 3056831 7359009 5486372 4840050 

482 

-7910400 -721410 -8912605 -3235810 -4996514 

-818 7i=12 

5456000 801567 5458450 4706456 4440571 *„= 5354228880 

909 

-2819200 -561097 -5820915 -294519 -6108400 

-636 

9127680 224438 8328366 117807 6443360 

254 

-7499520 -2008803 -5466306 -2444769 -3698464 

-44 

4841600 2905237 3881046 6773950 493200 

14642 

+.268946291E-007 

-192296598E-005 

+.448692063E-004 

-504778570E-003 

+.323058285E-002 

-128146453E-001 

+.329519451E-001 

-559006211E-001 

+.621118012E-001 

-434782609E-001 

+.173913043E-001 

-303030303E-002 

Jovial Program 

1.00 ARRAY A 16 17 I; 

2.00 ARRAY C 16 17 I; 

3.00 ARRAY V 17 10 I; 

4.00 ARRAY W 12 10 I; 

5.00 ARRAY MOD 10 I; 

6.00 ARRAY INV 10 I; 

7.00 ARRAY Z 10 I; 

8.00 ITEM SS F; ITEM TT F; ITEM NUM F; ITEM DEN F; 

9.00 VALUES MOD 9999889 9999901 9999907 9999929 9999931 9999937 

10.00 9999943 9999971 9999973 9999991; 

11.00 VALUES INV 7369185 2629157 6854848 9534445 6502552 

12.00 3472916 5202223 3251378 129200 5053746; 

13.00 VALUES W 4977739 8027834 5690123 538898 3121137 3768638 5584115 

14.00 115082 9999483 5278071 15291 9546516 9532061 5812518 

15.00 9202123 4699699 120098 9999471 6263753 7229845 9363502 8209581 

16.00 9936740 6614348 4222392 122714 9999465 7782099 9773154 8216183 

17.00 5103433 3253672 3499315 2247429 132922 9999443 3659841 431474 

18.00 2231991 6034197 1042630 5307831 2048927 133898 9999441 9722683 

19.00 2791024 2958205 5033108 173657 4212244 1433022 136874 9999435 

20.00 6535597 4607850 2424276 3272840 1038248 8705894 785509 139922 

21.00 9999429 2501001 296996 4616307 9611859 4673074 5857613 7307785 
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22.00 155097 9999401 9352927 6498696 8521281 1579588 4721737 3700720 

23.00 7028477 156241 9999399 8058781 9352296 9214244 9260206 9463057 

24.00 2740085 4312028 166897 9999381 7470971 3888109 6424092 5872382 

25.00 4092685 4802284 3931831 2809951 172468 9999372 8735485 1944054 

26.00 3212046 7936191 2046342 7401142 6965915 1404975 86234 4999686; 

27.00 READ M ;N=M + 1; 

28.00 PRINT 14H (MATRIX BY ROWS); 

29.00 1 = 0; 

30.00 Rl. H = 0; 

31.00 R2.READ C[I J]; 

32.00 IF J EQ N-2; GOTO R3; 

33.00 J = J + 1: GOTO R2; 

34.00 R3. IF I EQ M-l; GOTO R4; 

35.00 1 = 1 + 1; GOTO Rl; 

36.00 R4. PRINT 10H (RIGHT SIDE); 

37.00 1 = 0; 

38.00 R5. READ C[I ,N-1]; 

39.00 IF I EQ M-l; GOTO CI; 

40.00 1=1 + 1; GOTO R5; 

41.00 CI. I = 0;C2.J = 0;C3.V[I ,J] = 0;IF J EQ 9;GOTO C4; 

42.00 J = J + l;GOTO C3;C4.IF I EQ M ;GOTO R6;I = I + l;GOTO C2; 

43.00 R6. L = 0; 

44.00 R7. P=MOD[L]; 

45.00 1 = 0; 

46.00 R10. J = 0; 

47.00 R8. REMQUO (C[l J] ,P = Q ,A[I J]); 

48.00 IF I EQ N-l; GOTO R9; 

49.00 J = J + 1; GOTO R8; 

50.00 R9. IF I EQ M-l;GOTO Rll; 

51.00 I = I+l;GOTO R10; 

52.00 Rll. K = 0;SGN =1; 

53.00 Q0.R=K ;EX=1; 

54.00 Q1.B = A[R ,K] ;GOTO EO; 

55.00 Q2. IF D EQ l;GOTO Q3; 

56.00 IF R EQ M-l:GOTO ERR;R = R+ l;GOTO Ql; 

57.00 Q3.IF R EQ K ;GOTO Q5; 

58.00 T = K ;SGN = -SGN; 

59.00 Q4.F=A[R ,T] ;A[R ,T]=A[K ,T] ;A[K ,T] = F; 

60.00 IF T EQ M ;GOTO Q5; 

61.00 T = T + l;GOTO Q4; 

62.00 Q5.S = 0;Q8.IF K EQ S;GOTO Q7; 

63.00 G = X*A[S ,K] ;REMQUO (G ,P=Q ,G); 

64.00 T=K; 

65.00 Q6.REMQUO (G*A[K ,T] ,P = Q ,F);F=A[S ,T]-F ;REMQUO (F ,P = Q ,A[S ,T]) ; 

66.00 IF T EQ M ;GOTO Q7;T = T+ l;GOTO Q6; 

67.00 Q7. IF S EQ M-l ;GOTO Q9;S = S + 1 ;GOTO Q8; 

68.00 Q9. IF K EQ M-l;GOTO Q12;K=K+l;GOTO Q0; 
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69.00 Q12.K = 0;EX = 2;Q15.B = A[K ,K] ;GOTO E0; 

70.00 Q13.F = X *A[K ,M] ;REMQUO (F ,P = Q ,F); 

71.00 IF F GQ 0;GOTO Q14;F = F + P; 

72.00 Q14. A[K ,M] = F; 

73.00 IF K EQ M-1;G0T0 P16;K=K+ 1;G0T0 Q15; 

74.00 Q16.G0T0 L18; 

75.00 L18.I = 0;R12.SGN=SGN *A[I ,1] ;REMQUO (SGN ,P=Q ,SGN); 

76.00 IF I EQ M-l;GOTO F3;I = I+ l;GOTO R12;F3. IF SGN GQ 0; 

77.00 GOTO F1;SGN = SGN + P;GOTO F1;F1.K = M-1;F2.T = A[K ,M] *SGN; 

78.00 REMQUO (T ,P = Q ,A[K ,M]);IF K EQ 0;GOTO R13;K = K-l;GOTO F2; 

79.00 R13.EX=1;K = 0;REMQUO (INV[L] *SGN ,P = Q ,R); 

80.00 R14.Z[K] = W[L ,K] *R + V[M ,K] ;IF K EQ 9;GOTO NRM ;K=K+1;GOTO R14; 

81.00 R15.K = 0;R16.V[M ,K] = Z[K] ;IF K EQ 9;GOTO R17;K = K+ l;GOTO R16; 

82.00 R17.EX = 2;I = 0;R20.K=0;REMQUO (INV[L] *A[I ,M] ,P=Q ,R); 

83.00 R18.Z[K]=W[L ,K] *R + V[I ,K] ;IF K EQ 9;GOTO NRM ;K=K+l;GOTO R18; 

84.00 R22.K=0;R19.V[I ,K] = Z[K] ;IF K EQ 9;GOTO R21;K=K+ l;GOTO R19; 

85.00 R21.IF I EQ M-l;GOTO RM ;I = I + l;GOTO R20; 

86.00 RM. IF L EQ 9;GOTO RR ;L = L+l;GOTO R7; 

87.00 RR.I = 0;RR 1 .K=9;CMP .IF V[I ,K] LS W[ll ,K] ;GOTO RR2; 

88.00 IF V[I ,K] GR W[ll ,K] ;GOTO RR3;IF K EQ 0;GOTO RR3; 

89.00 K=K-l;GOTO CMP; 

90.00 RR3.EX = 3;K = 0; RR4.Z[K] = V[I ,K]-W[10 ,K] ;IF K EQ 9;GOTO NRM; 

91.00 K=K+l;GOTO RR4; 

92.00 RR6.K = 0;RR5.V[I ,K] = Z[K] ;IF K EQ 9;GOTO RRl;K = K+l;GOTO RR5; 

93.00 RR2.IF I EQ M ;GOTO PR ;I=I+l;GOTO RR1; 

94.00 PR.I = 0; 

95.00 PR1.PRINT V[I ,0] ,V[I ,1] ,V[I ,2] ,V[I ,3] ,V[I ,4]; 

96.00 PRINT V[I ,5] ,V[I ,6] ,V[I ,7] ,V[I ,8] ,V[I ,9]; 

97.00 IF I EQ M ;GOTO G4;I = I+1 ;GOTO PR1; 

98.00 NRM.K = 0; 

99.00 Nl.REMQUO (Z[K] ,10000000 = Q ,Z[K]) ;Z[K+ 1] = Q + Z[K+1]; 

100.00 IF K EQ 8;GOTO N2;K=K+l;GOTO Nl; 

101.00 N2.K = 9; 

102.00 N3.IF Z[K] GR 0;GOTO N6;IF Z[K] LS 0;GOTO N7; 

103.00 IF K EQ 0;GOTO NN ;K=K-l;GOTO N3; 

104.00 N6.K = 0;N4.IF Z[K] GQ ; GOTO N5;Z[K] = Z[K] + 10000000; 

105.00 Z[K+1]=Z[K+1]-1;N5.IF K EQ 8;GOTO NN;K=K+ l;GOTO N4; 

106.00 N7.K =0;N8.IF Z[K] LQ 0;GOTO N9;Z[K] = Z[K]- 10000000; 

107.00 Z[K+1]=Z[K+1]+1 ;N9.IF K EQ 8;GOTO NN ;K=K+l;GOTO N8; 

108.00 NN.IF EX EQ l;GOTO R15;IF EX EQ 2;GOTO R22;IF EX EQ 3; 

109.00 GOTO RR6; 

110.00 E0.F = B ;D = P ;S=1 ;X = 0; 

111.00 El.REMQUO (F ,D = G ,T); 

112.00 IF T EQ 0;GOTO E2; 

113.00 F=D ;D = T ;T=S-G *X; 
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114.00 S = X ;X = T ;GOTO El; 

115.00 E2.IF D GQ 0;GOTO E3; 

116.00 D = -D ;X = -X; 

117.00 E3.IF X LS 0;X = P + X ;IF EX EQ l;GOTO Q2;IF EX EQ 2;GOTO Q13; 

118.00 ERR .PRINT P ,13H (DIVIDES DET A) ;GOTO RM; 

119.00 G4.I = 0;G3.NUM = V[I ,9] ;DEN = V[M ,9] ;K = 2; 

120.00 G1.TT = V[I ,10 -K] ;SS = V[M ,10 -K] ;NUM= 10000000. *NUM + TT; 

121.00 DEN=10000000. *DEN+SS ;IF K EQ 10;GOTO G2;K=K+1 ;GOTO Gl; 

122.00 G2.PRINT NUM/DEN ;IF I EQ M-l;STOP ;I = I+1 ;GOTO G3; 



(Paper 71B4-240) 
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